
/**
  ******************************************************************************
  * Copyright 2021 The grapilot Authors. All Rights Reserved.
  * 
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.
  * You may obtain a copy of the License at
  * 
  * http://www.apache.org/licenses/LICENSE-2.0
  * 
  * Unless required by applicable law or agreed to in writing, software
  * distributed under the License is distributed on an "AS IS" BASIS,
  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  * See the License for the specific language governing permissions and
  * limitations under the License.
  * 
  * @file       Atmtable.c
  * @author     baiyang
  * @date       2021-7-14
  ******************************************************************************
  */

/*----------------------------------include-----------------------------------*/
#include <stdio.h>
#include <string.h>
#include <math.h>

/*-----------------------------------macro------------------------------------*/
#define TZERO    288.15f                  /* sea-level temperature, kelvins */
#define SUTH     110.4f                  /* Sutherland's constant, kelvins */
#define BETAVISC 1.458E-6f                           /* viscosity constant */

/*----------------------------------typedef-----------------------------------*/

/*---------------------------------prototype----------------------------------*/

/*----------------------------------variable----------------------------------*/

/*-------------------------------------os-------------------------------------*/

/*----------------------------------function----------------------------------*/
void Atmosphereconst(float  alt,              /* geometric altitude, km.  */
    float *sigma,                                /* density/sea-level standard density */
    float *delta,                            /* pressure/sea-level standard pressure */
    float *theta)
{
    const float GMR = 34.163195f;    //理想气体常数
    const float REARTH = 6356.766f;

    int i, j, k;
    float h, tgrad, deltah, tbase, tlocal;

#define NTAB 8
    static float htab[NTAB] = { 0.0f, 11.0f, 20.0f, 32.0f, 47.0f,
        51.0f, 71.0f, 84.852f };
    static float ttab[NTAB] = { 288.15f, 216.65f, 216.65f, 228.65f, 270.65f,
        270.65f, 214.65f, 186.946f };
    static float ptab[NTAB] = { 1.0f, 2.2336110E-1f, 5.4032950E-2f, 8.5666784E-3f,
        1.0945601E-3f, 6.6063531E-4f, 3.9046834E-5f, 3.68501E-6f };
    static float gtab[NTAB] = { -6.5f, 0.0f, 1.0f, 2.8f, 0.0f, -2.8f, -2.0f, 0.0f };

    h = alt*REARTH / (alt + REARTH);        /* geometric to geopotential altitude  */

    i = 0; j = NTAB - 1;                      /* starting values for binary search */
    do
    {
        k = (i + j) / 2;
        if (h < htab[k]) j = k; else i = k;
    } while (j > i + 1);

    tgrad = gtab[i];
    tbase = ttab[i];
    deltah = h - htab[i];
    tlocal = tbase + tgrad*deltah;
    *theta = tlocal / ttab[0];                            /*  temperature ratio  */

    if (tgrad == 0.0f)                                    /*  compute pressure ratio  */
        *delta = ptab[i] * expf(-GMR*deltah / tbase);
    else
        *delta = ptab[i] * powf(tbase / tlocal, GMR / tgrad);

    *sigma = *delta / *theta;
}

void SimpleAtmosphere(
    const float  alt,                      /* geometric altitude, km.  */
    float *sigma,                /* density/sea-level standard density */
    float *delta,              /* pressure/sea-level standard pressure */
    float *theta)        /* temperature/sea-level standard temperature */

                         /* Compute the temperature,density,and pressure in the standard atmosphere */
                         /* Correct to 20 km.  Only approximate thereafter.    */
{
    float REARTH = 6356.766f;               /* polar radius of the Earth (km) */
    float GMR = 34.163195f;                                /* gas constant */

    float h;

    h = alt*REARTH / (alt + REARTH);        /* geometric to geopotential altitude  */

    if (h<11.0f)
    {                                                       /* Troposphere */
        *theta = (288.15f - 6.5f*h) / 288.15f;
        *delta = powf(*theta, GMR / 6.5f);
    }
    else
    {                                                      /* Stratosphere */
        *theta = 216.65f / 288.15f;
        *delta = 0.2233611f*expf(-GMR*(h - 11.0f) / 216.65f);
    }

    *sigma = *delta / *theta;
}   /* ---------------------------------- End of function SimpleAtmosphere */

    /* ======================================================================= */
float MetricViscosity(const float theta)
{
    float t = theta*TZERO;                           /* temperature in kelvins */
    return BETAVISC*sqrtf(t*t*t) / (t + SUTH);         /* viscosity in kg/(m-sec) */
}   /* ----------------------------------- End of function MetricViscosity */

/*------------------------------------test------------------------------------*/


